Atmosphere‐Snow Exchange Explains Surface Snow Isotope Variability

Abstract The climate signal imprinted in the snow isotopic composition allows to infer past climate variability from ice core stable water isotope records. The concurrent evolution of vapor and surface snow isotopic composition between precipitation events indicates that post‐depositional atmosphere‐snow humidity exchange influences the snow and hence the ice core isotope signal. To date, however, this is not accounted for in paeleoclimate reconstructions from isotope records. Here we show that vapor‐snow exchange explains 36% of the summertime day‐to‐day δ18O variability of the surface snow between precipitation events, and 53% of the δD variability. Through observations from the Greenland Ice Sheet and accompanying modeling we demonstrate that vapor‐snow exchange introduces a warm bias on the summertime snow isotope value relevant for ice core records. In case of long‐term variability in atmosphere‐snow exchange the relevance for the ice core signal is also variable and thus paleoclimate reconstructions from isotope records should be revisited.


Introduction
This Supplementary Material includes some detailed complementary information to the main text that may be of interest to some readers. The two field campaigns are compared in terms of vapor and snow isotopic composition in Table S1 and Figure S1. The calibration protocol for the CRDS water vapor isotope measurements is given in Text 1. Details on the Eddy-Covariance systems are given in Text 2. The noise level estimates used in the Monte Carlo simulations for the model uncertainty estimate are outlined in Text 3. Text 4 and Figure S5 explain how we estimated the uncertainty of the snow surface isotopic composition. The snow surface model is described in more detail in Text 5 and Figure S2 while the results of the model sensitivity study are shown in Figure S4 and detailed values are given in Table S2. Finally, the results of the experiment using equilibrium fractionation (ii) during the sublimation process are shown in Figure S3.

Text 1: Calibration protocol for CRDS measurements of water vapor
The CRDS measurements were corrected for the instrument's humidity dependency and calibrated against the VSMOW-SLAP scale following the approach of Steen-Larsen et al. (2013). Calibrations were performed every day using vaporized secondary standards at a constant humidity level. Secondary standards were provided by the Institute of Arctic and Alpine Research, University of Colorado and bracketed the observed range of dvalues. Since the CRDS was observed to be stable within the individual field campaigns, all calibration runs were grouped together to obtain one set of VSMOW-SLAP calibration curves for dD and d 18 O per campaign. The CRDS' humidity dependency was corrected using humidity calibrations carried out several times throughout the campaigns. The humidity dependency curves were established for each standard by providing constant vapor streams at different humidity levels. The humidity dependency response curves for each campaign were linearly interpolated in the d-space to account for the dependency of the response curves on the standard's isotopic composition. The humidity mixing ratio of the CRDS was calibrated against a custom-made dew-point generator for a temperature range between -40°C and 0°C in the laboratory.

Text 2: Specifications of the eddy-covariance systems
In 2018, the system consisted of a CSAT3 sonic anemometer and a Krypton Hygrometer (KH20) both from Campbell Scientific (CS). Latent heat flux (LE) estimates for 2018 were calculated for averages of 10min using a custom processing software and postprocessed to 30min means. In 2019, 30min resolution surface fluxes were measured using a CS Irgason system and the processing software EasyFlux which follows established correction protocols (e.g. Mauder et al., 2006).

Text 3: Uncertainty estimate of modeled d using Monte Carlo simulations
For estimating the uncertainty of the modeled d*, we performed 1000 Monte Carlo simulations for each day-to-day modeling period by adding random noise to the input parameters LE, dV, initial dS and atmospheric humidity q. For LE, we account for EC measurement accuracy errors by multiplying the input timeseries with a factor drawn from a normal distribution around a mean of 1 and SD 20%. The parameters h and Ts are calculated from LE subsequently. For the other input parameters, we account for the uncertainty of the measurements by adding noise to the individual datapoints drawn from a normal distribution with mean 0 and a representative SD for the respective parameter, i.e. SD of dv is 0.23‰ for d 18 O and 1.4‰ for dD (Steen-Larsen et al., 2013). Due to the absence of a manufacturer precision statement for the CRDS humidity readings, we infer a relative error of 0.2% based on the calibration pulses. Finally, we define the uncertainty of the estimated mean snow surface isotopic composition dS from an additional snow sample set as 0.53‰ and 3.8‰ for d 18 O and dD, respectively (Text 1). Since d* is the subtraction of two consolidated snow samples we propagate the error to find the uncertainty on the observed d* to be 0.74‰ for d 18 O and 5.4‰ for dD which are indicated as error bars in Figure 3 in the main text.
Text 4: Estimating the uncertainty of the average 0.5cm snow-layer isotopic composition To estimate the uncertainty associated with the consolidated snow surface isotopic composition sample as representative value for the mean snow surface isotopic composition we make use of snow samples collected from ten positions along a 90m long transect perpendicular to the prevailing wind direction. Snow of the upper 0-0.5cm layer was sampled with a 10m distance between sampling locations. The sampling scheme was repeated on 30 days in 2019. We calculate the deviation of each individual sample from the respective sample mean and calculate the standard deviation of all data points. The distribution is shown as histogram in Fig. S5. This yields an uncertainty on the consolidated snow isotopic composition of e=0.53‰ for d 18 O and e=3.8‰ for dD.

Text 5: Snow surface model
The evolution of the snow surface isotopic composition is simulated by setting up a snow surface model using a mass balance approach. At initializing, a H=5 cm thick snowpack domain is divided into ten layers with j = 1 being the surface layer in contact with the atmosphere. The mass (m) per unit surface area of the surface layer is calculated using the daily observed snow densities (ρ ! ).

Vapor-Snow Exchange
At every time step (i), with dt=30 min resolution, an amount of mass (dm) is either removed or added to the snowpack during sublimation (LE > 0) or deposition (LE < 0) conditions, respectively. The LE measurements together with the latent heat of sublimation (L " = 2.834e # J kg $% ) are used to calculate the sublimation rate (s in kg m $& s $% ), from which dm is calculated following: (2) The mass balance per time step equals: With R ! and R ) being isotopic ratio of snow and humidity flux, respectively.
From this we can find the partial time (t) derivative for small mass changes: Converted to δ-notation in ‰ this equals: The isotopic composition of the sublimation flux is computed in analogy to the evaporation flux as formulated by Merlivat and Jouzel (1979): Where h is the saturation deficit, i.e. (1 − h) is the humidity (q) gradient inducing the humidity flux. δ , and q , are the isotopic composition of the vapor and the humidity at the 2 m level. The surface is assumed to be at saturation q "-. (T ! ) with regards to the surface temperature T ! . T ! , and hence q "-. and h, are calculated from LE and q , solving a humidity flux bulk equation following Van As (2011) with stability correction functions for stable (Holtslag & De Bruin, 1988) and for unstable (Paulson, 1970) conditions. We use the median of observed roughness length (z0) values during neutral atmospheric stability in 2019 as constant momentum roughness length of z0=1.3·10 -4 m and calculate the scalar roughness length from it following Andreas (1987). α *+ (Ts) and k are the equilibrium and kinetic fractionation coefficients, respectively with k !" / = 6 ‰ and k 0 = 0.88 ⋅ k !" / . Icevapor equilibrium fractionation factors of Majoube (1971) and Merlivat & Nief (1967) are used for the sublimation flux and any other occasion of equilibrium fractionation. During times of 0.98 < h < 1 the snow isotopic composition is not changed by the sublimation flux since δ ) approaches unrealistic values.
During deposition conditions, we use the empirically derived linear functions of Wahl et al. (2021) to calculate the humidity flux isotopic composition from the 2m vapor isotopic composition.: The vapor-snow humidity exchange only affects the isotopic composition of the surface layer and its mass. Additionally, a diffusion model is realized, that modulates the isotopic composition along the gradients within the air space of all layers within the 5 cm snowpack at every time step.

Isotope Diffusion in Snow
Following Johnsen et al. (2000), the diffusivity (Ω) of water vapor in air (a) in cm & s $% is given by reference (P 1 = 1 atm) and ambient pressure (P) and absolute (T 1 = 273.15 K) and snow temperature (T). T of the entire 5cm snowpack model domain is calculated for every time step as average from in-situ observed snow temperatures at 10cm depth and calculated Ts.
Ratios of isotope diffusivities are taken from Merlivat (1978) which yields diffusivities of the isotope species (Ω - * ) in air with: In snow, the pore space shape limits the diffusion, which can be described with the tortuosity factor (τ), that is assumed to be a function of density (Schwander et al., 1988): with density of ice ρ ' = 917 kg m $7 . The diffusivity in firn (here snow) can then be calculated with molar mass of water (mo), saturation vapor pressure over ice (p) from Goff & Gratch (1946) and the gas constant R = 8.314 Jmol $% K $% : Diffusion of isotopes within the snow pack can therefore be described by:

Variable Snow Height
To account for mass transfer and resulting changes in snow pack height, we introduce the dimensionless height (ξ), defined by the depth profile (z) and H: This way the mass gain or loss are equally distributed over all snowpack layers which preserves a surface layer close to the original 0.5cm and thus prevents an artificial enhancement of the sublimation effect in case of high net sublimation, i.e. high mass removal which would otherwise only thin the top layer.

Full Snow Surface Model
The and B for surface flux influences on top layer j=1: = 1 1 0 ⋮ 0 2 then reads for either isotope species δ * and timestep i:

Fig. S1 The surface and sub-surface snow isotope variability
Timeseries of snow isotopic composition in different depths are presented for the 2018 and 2019 field season in (a) and (b), respectively. Note that the values presented are calculated from the snow samples (available on PANGAEA) that were sampled from the surface to the sampling depth under the assumption of homogenous mixing of the different layers.

Fig. S2 Scheme of the snow surface model
The snowpack model domain is shown where the surface layer (j1) is influenced by the humidity flux (LE) and diffusion between snow layers (j). Model input parameters that are updated from observations are marked with *, while the other parameters are calculated from these observations. The snowpack isotope profile is initialized at the time of a snow sampling event and evolves driven by the humidity flux and diffusion until the time of the next sampling event the day after. The modeled change ( d*) in the top snow layer's isotopic composition is then compared to the observed change.   . The x-axis shows the SD of the distribution of d* resulting when perturbing the input parameters of the model as shown in Table S2. The y-axis shows the modeled change in surface isotopic composition using the default input parameters. The median of all SD per disturbed parameter are shown on the lower axes and given in Table S2.